Thermal collapse of a granular gas under gravity 
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Free cooling of a gas of inelastically colliding hard spheres represents a central paradigm of kinetic 
theory of granular gases. At zero gravity the temperature of a freely cooling homogeneous granular 
gas follows a power law in time. How does gravity, which brings inhomogeneity, affect the cooling? 
We combine molecular dynamics simulations, a numerical solution of hydrodynamic equations and 
an analytic theory to show that a granular gas cooling under gravity undergoes thermal collapse: it 
cools down to zero temperature and condenses on the bottom of the container in a finite time. 

PACS numbers: 45.70.Qj, 47.70.Nd 



Granular a' low-density fluid of inelastic hard 

spheres, is a simple model of granular flow, and it has 
attracted much attention from physicists Q, Q . An un- 
driven granular gas loses its kinetic energy via inelastic 
collisions. In the Homogeneous Cooling State (HCS) the 
temperature T of a dilute granular gas decays according 
to Haff's law [![, T{t) = T (l + t/t )- 2 , where in two di- 
mensions to = \Jit/2 (1 — r 2 ) dnoT^ is the cooling time, 
rip is the (constant) number density of the particles, d is 
the particle diameter and r is the coefficient of normal 
restitution. The HCS, and deviations from it, provide a 
rich testing ground for the ideas and methods of kinetic 
theory of granular gases, and it has been investigated 
in many theoretical works, see Ref. Q and references 
therein. Direct experimental observation on the HCS is 
difficult, not the least because of gravity. Therefore it is 
somewhat surprising that there have been no theoretical 
studies of the effect of gravity on the free cooling of a 
granular gas. It is intuitively clear that gravity forces 
grains to sink to the bottom of the container, where in- 
creased density enhances the collision rate and causes 
"freezing" of the granulate. However, no quantitative 
analysis of this process has ever been performed. Here 
we combine molecular dynamics (MD) simulations, a nu- 
merical solution of granular hydrodynamic equations and 
analytical theory to develop a detailed quantitative un- 
derstanding of this cooling process. Our main result is 
that, in a striking contrast to Haff's law, the gas under- 
goes thermal collapse: it cools down to zero temperature 
and condenses on the bottom plate in a finite time ex- 
hibiting, close to collapse, a previously unknown univer- 
sal scaling behavior. 

MD simulations. We employed an event-driven algo- 
rithm Q to simulate a free cooling of an initially dilute 
gas of N ^> 1 identical nearly elastic, 1 — r <C 1, hard 
disks of unit diameter and mass in a two-dimensional 
container of width L x and infinite height. The (elastic) 
bottom of the container is at y = 0, the (elastic) side 
walls are at x = and L x . L x is chosen small enough 
so that any macroscopic structure in the lateral direc- 
tion is suppressed. The gravity acceleration g acts in the 
negative y direction. Figure shows four snapshots of 



FIG. 1: Snapshots of an event-driven MD simulation at indi- 
cated times for N = 5642, L x = 10 2 , r = 0.995, T = 10 and 
g = 0.01. Only a part of the box is shown. 

a typical simulation where, at t = 0, particles have a 
Maxwell velocity distribution, and a Boltzmann density 
profile at constant temperature To Q. Collapse of all 
particles to the bottom is observed at time t c = 7770. 
The circles in Fig. show the time history of the sim- 
ulated total kinetic energy of the gas, normalized to its 
value at t — 0. One can see that the total energy drops 
to zero in a finite time. We observed a similar behavior 
in a wide range of parameters, and also for a different, 
non-isothermal initial state, prepared by replacing the 
elastic bottom plate by a "thermal" bottom plate [3] and 
waiting until a steady state is reached. In the latter case 
the initial transient is somewhat different, but the energy 
decay law close to collapse remains the same, see Fig. 

Hydrodynamic theory. The observed energy decay dy- 
namics are remarkably captured by hydrodynamic equa- 
tions for the number density n(y,t), vertical velocity 
v(y,t) and granular temperature T(y,t). These equa- 
tions are systematically derivable from the Boltzmann 
equation generalized to account for inelastic collisions of 
hard disks 0. We assume a dilute gas, an assumption 
which becomes invalid close to collapse. Following Ref. 
, we rescale the variables using the gravity length scale 
A = To/g and the heat diffusion time td, = e^ 1 {\/ g) 1 / 2 . 
The scaled parameter e = ir~ 1 / 2 (L x /Nd) <C I is of the 
order of the inverse number of layers of grains which form 
after the particles settle on the bottom. The smallness 
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FIG. 2: (a) Total kinetic energy of the grains, normalized 
to its value at t = 0, vs. time for the run in Fig. Q (cir- 
cles). Also shown are the results from the full hydrodynamic 
model Q-Q (black solid line), from the oj-equation j^J (red 
dashed line) , and from a MD simulation with a different, non- 
isothermal initial condition (black dash-dotted line). The re- 
spective hydrodynamic parameters are e = 10" 2 and A = 5. 
The inset shows the square root of the energy close to thermal 
collapse, (b) a space-time plot of T from the full hydrody- 
namic model for the same simulation. 

of e guarantees that td is much longer than the fast hy- 
drodynamic time tf — (X/g) 1 ^ 2 . We measure v in units 
of vq = X/td, n in units of uq = N/XL X , and T in units 
of To. Furthermore, we exploit the one-dimensionality 
of the flow and go over to Lagrangian mass coordinate 
to = L n{y' ,t)dy' which varies between at the bottom 
and 1 (the total rescaled mass of the gas) as y — > oo. The 
resulting rescaled hydrodynamic equations are @ : 

d t (l/n) = d m v, (1) 
s 2 d t v = -d m {nT) - 1 + (e 2 /2) d m (nT^ 2 d m v) , (2) 
d t T + nTd m v = (e 2 /2) nT x l 2 {d m vf + 
(4/3) d m (nd m T 3 / 2 ) - 4A 2 toT 3 / 2 , (3) 

In addition to e, Eqs. ©-© include the parameter A 2 = 
■ a which shows the relative role of the inelastic energy 
loss and heat diffusion. At the boundaries y — and oo 
we demand zero fluxes of mass, momentum and energy 
0, which yield v — d m T = at m = and nd m v = 
nd rn T — at to = 1 . 

We solved Eqs. OJ~0 numerically in a wide range 
of parameters, using a variable mesh/variable time step 
solver . The blue solid line in Fig. |2t depicts the sum 
of the thermal energy and macroscopic kinetic energy of 
the gas versus time for the same parameters and initial 
condition as in the MD simulation indicated by the cir- 
cles. The agreement is excellent, and thermal collapse is 
clearly observed 6] . At the very early stage of the cooling 
[with duration of 0(r/)] we observed shock waves which 
form at large heights, cause a transient heating of the gas 
there, and escape to m = 1 (y = oo), see Fig. [2fc>. 

Quasi-static flow. If e <C min(l,A~ 2 ) then, after the 
brief transient, a quasi-static flow sets in. Here the e 2 - 
terms in Eqs. @ and can be neglected, and Eq. 
@ reduces to the hydrostatic condition d m {nT) + 1 = 
which yields nT = 1— to. Substituting n = (1— m)/T into 
Eq. © and using Eq. (JIJ, we obtain a closed nonlinear 
equation for a new variable u>(m,t) — T^'^{m,€): 

(ddtLo — d m [(1 — m)d m u>] — A 2 (l — m)uj , (4) 
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FIG. 3: Hydrodynamic fields predicted by the u-equation for 
A = 1: T(m,t) at times separated by At = 0.2 (a), and T(y,t) 
(b), v(y,t) (c), and T(m,t)(t c — t)~ 2 (d) at the same times. 
The late-time curves in d collapse into a single curve. 

We will call Eq. J3J the w-equation; it was derived, in 
another context, in Ref. 0. 

We solved the w-equation numerically [with the no- 
flux boundary conditions (1 — m)d m u> = at m = 
and 1], using the same solver 0. A typical example is 
shown in Fig. [5^,. Here we launched the computation 
at scaled time t — 0.04 when the hydrostatic condition 
nT = l — m already holds well, and used the temperature 
profile, computed with the full hydrodynamic solver, as 
the initial condition. One can see that the w-equation 
provides a faithful description of the later stage of the 
cooling. Figure|3]shows a different example of the cooling 
dynamics, as described by the w-equation starting from 
w(m, t = 0) = 1. Here we show the T- and i>-profiles in 
both Lagrangian and Eulcrian coordinates. In all simu- 
lations thermal collapse is observed at a time t c which 
goes down as A increases. The collapse occurs simulta- 
neously on the whole Lagrangian interval (0, 1), see Fig. 
[3^,. As the density n = (1 — m)/T blows up at t = t c at 
all to £ [0, 1), this Lagrangian interval corresponds to a 
single Eulerian point y — 0. Therefore, at time t = t c all 
of the gas condenses on the bottom plate y = and cools 
to a zero temperature [To| . 

Separable solution close to collapse. As Fig. [3Jl im- 
plies, uj(m,t) becomes separable as t — * t c . This remark- 
able solution can be written as 

cj(m,t) = (t e -t)Q(m), (5) 

where Q(m) is determined by the nonlinear ODE 

[(1-to)0']'-A 2 (1-to)0 + Q 2 = (6) 

(the primes denote m-derivatives) and the boundary con- 
ditions (1 — to) Q' — at to = and 1. Function Q(m) is 
uniquely determined by A and, at fixed A, can be found 
numerically by shooting. In addition, we found Q(m) 
perturbatively for small and large A: 
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I. A 2 <C 1. As it can be checked a posteriori, in this 
case Q(m) = 0(A 2 ). Furthermore, as the heat diffu- 
sion dominates over the inelastic energy loss, the so- 
lution must be almost constant on the whole interval 
< m < 1. Therefore, we seek a solution in the form 
Q(m) = A 2 Qo + A 4 <2i(m) + A 6 Q 2 (m) + . . . . Substituting 
this in Eq. 10 and equating terms of the same order in 
A 2 , we obtain the asymptotic solution 



Q(m) 



A 2 . 
2 



A 4 
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C(A 6 



(J) 



We checked that this solution is in excellent agreement 
with numerical solutions of the w-equation at small A. 

II. A 2 ^> 1. Here it is convenient to stretch the La- 
grangian coordinate, £ = A(l — to), and time r = At, so 
that A drops from the w-equation 



u)d T ui = £<9|w + d^uj — £a 



(8) 



but enters the integration interval (0, A], whereas the 
boundary condition are £,d^oj — 0. The separable so- 
lution is c<j(£, t) = (t c — T)g(£), while the boundary-value 
problem for g(£) is 



(&')' - £<z + 1 2 = o, £<z' = o at £ = 0,A. 



(9) 



At £ 3> 1 g(£) is exponentially small, so one can 
drop the q 2 -term and obtain ?&(£) = C[ATo(£) + 
i^ 1 (A)/r 1 (A)/o(e)] (where K 0)1 (£) and / ,i(0 are the 
modified Bessel functions), which obeys the boundary 
condition at £ = A. This solution with C = 0.951 agrees 
well with the full numerical solution already at £ > 1, 
see Fig. and therefore is valid everywhere except the 
thin boundary layer at m — ► 1. As the q 2 -term originates 
from the u>d T ui term in Eq. (JSJ, we realize that, almost 
everywhere, the energy loss at late times is balanced by 
the heat conduction, while the boundary layer at m — > 1 
serves as a dynamic "bottleneck" of the cooling. 

Outside of a thin boundary layer near £ = A (or 
to = 0), the solution is close to the solution of the same 
equation but on the semi- infinite interval (0, oo). The 
latter one, 9oo(£), is parameter-free and can be found nu- 
merically. The shooting starts at the left boundary £ = 
which is a regular singular point of Eq. JjJJl . We demand 
that <?£,(£ = 0) be finite, which yields g^(0) = -goo(0) 2 . 
The shooting procedure gives a unique value of <7oo(0) for 
which the solution does not diverge toward +oo or — oo 
at large £. We find qo = <7oo(0) = 1.633356 . . . ; the re- 
spective asymptotic profile (?oo(£) is the envelope of the 
numerical profiles g(£) for different A in Fig. HJ,. 

Early dynamics and collapse time. To find the collapse 
time t c [a free parameter in the separable solution JjjJ], 
one needs to solve the w-equation with a given initial 
condition. In general, this can only be done numerically. 
We obtained analytic estimates of t c , separately for small 
and large A, for uj{m, t = 0) = 1. 

For A< 1 the separable solution (JSJ and Q is valid, at 
order A 2 , at any t > 0. This yields the leading-order esti- 
mate t c ~ 2A~ 2 . For A> 1 the initial stage of the cool- 
ing dynamics should be addressed separately. We notice 
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FIG. 4: (a) Numerical solutions of Eq. © (symbols) and the 
corresponding "bulk solutions" <#,(£) (lines) for A = 5, 10, 15. 
(b) The implicit solution 1|'I0I (the blue lines) versus the nu- 
merical solution of the tj-equation (the black lines) for A = 8.0 
at three Lagrangian points: m = 0,0.5 and 1.0. The red 
straight line shows the smoothly matched asymptotic separa- 
ble solution at m = 1 (b). 

that at early times the term £<9 2 w in Eq. (JSJ) is small com- 
pared to the rest of terms. With this term neglected Eq. 
© reduces to a first-order equation, uid T uj — d^uj = — 
which is soluble by characteristics. The solution a>(£, r), 
in an implicit form, is 



V2e-« 2 /M£,r) 



I 



dz = 



(10) 



it is depicted, at points m = 0,1/2 and 1, in Fig. 0Jd. 
At £ 3> 1 (where most of the gas is located), Eq. (fT0|) 
predicts an early-time asymptote w(£,r -C 1) = 1 — £t 
which can be also obtained directly from Eq. (JSJ) with the 
heat conduction neglected completely. The "bottleneck" 
of cooling, however, is at large heights, £ <C 1, where 
the gas is very dilute. An early-time asymptote there, as 
predicted from Eq. I|10|l . is 

uj ~ 1 - £r - t 2 /2 = 1 - A 2 (l - m)t - A 2 t 2 /2 . (11) 

The implicit solution llOfl breaks down, at a given £, at 
time r ~ min (1, £ _1 ), and then the full w-equation must 
be solved. Eventually, as r approaches r c , the separable 
solution (J5J emerges. 

We stress that, at A ^> 1, the cooling process is highly 
nonuniform, see Fig. [SJi and b. For example, at m = a 
rapid initial decay co = 1 — A 2 i crosses over, after a short 
time t ~ A -2 , into a very slow decay u> = Aq{A)(t c — t), 
as q{A) is exponentially small. Meanwhile, at m = 1 a 
slow initial decay co = 1 — r 2 /2 crosses over into a rapid 
decay u> = qo(r c — r). At t = t c u> vanishes at all £ (that 
is, at y = 0). Note that the dynamics at to = 1 are 
independent of A in the stretched time t = At 

A good estimate of the collapse time r c at large A 
can be obtained by matching, at £ = 0, the late-time 
asymptote (r c — r)qo with the early-time solution i|10fl . 
see Fig. 0Jd. This yields an algebraic equation for r c : 



( 7 r/2) 1 / 2 erfi(2- 1 / 2r -i : 



= r c exp[l/(2r 2 )] (12) 



which has a unique solution t c ~ 1.10, or t c ~ 1.10A . 
For comparison, a numerical solution of Eq. (jSJ, for 
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large A, with the initial condition t = 0) = 1 yields 
t c = 1.145 which agrees with our estimate to about 4%. 
The slope of the collapsed curves w(£ = 0, r) for large A 
(Fig. 2]d) near r = r c , <9 T o;(£ — 0,t c ) = 1.624, is in very 
good agreement with the asymptotic value go = 1.633356. 
Our predictions of the A-dependence of t c are summa- 
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FIG. 5: (a) u(m = 0, t) at A = 10. The inset shows a blow-up 
near t c , and the asymptote (r c — r)g(A) (the red dashed line), 
(b) uj(m — l,t) vs. r = At for A = 5, 10, 15 collapse into a 
universal curve. 
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FIG. 6: The collapse time versus A. Numerical results from 
the adequation and from the full hydrodynamic equations Q- 
© with e = 0.025 are shown by symbols, the small- and 
large- A asymptotes are denoted by lines. 

rized in Fig. The small- and large-A asymptotes are in 
excellent agreement with numerical results. Returning to 
the dimensional units, we observe that, at A <C 1 the col- 
lapse time t c is much longer than the heat diffusion time. 
At A > 1 t c is of the order of (eA)" 1 ^ - (1 - r 2 )- 1 / 2 t f 



which, for nearly elastic collisions, is much longer than 
the free fall time t / . 

Having found r) we can find the rest of hy- 
drodynamic fields. Here we present the results for 
A 3> 1. In the early stage of cooling the gas den- 
sity is ti(£,t) = (£/A)(l — £t)~ 2 and velocity v(£, r) = 
A(A - £)[(A + £)r - 2]. Close to collapse the density 
n(Z> T ) = t( T c - t)- 2 <T 2 (£) blows up as (r c - r)" 2 . The 
gas velocity is v(£, r) = -2A(r c - r) / f V (£')/£'] d?. At 
r < t c it diverges logarithmically at £ = (that is, lin- 
early at y — > oo), but vanishes everywhere at r = r c , 
while the mass flux nv blows up at t = r c . Going back 

to Eulcrian coordinate y = (t c — t) 2 J^[q 2 (C)/C] d£', we 
see that the velocity field is simply v(y, t) — —2yj (t c — t). 

Summary. Our MD simulations and hydrodynamic the- 
ory depict a coherent picture of thermal collapse which 
develops in the process of a free cooling of a granular 
gas under gravity. One of the signatures of this pic- 
ture is the universal scaling behavior of the total energy 
E(t) ~ (t c - t) 2 as t -> t c . 

It would be interesting to test the quantitative predic- 
tions of our theory in experiment. A possible experiment 
can employ metallic spheres rolling on a slightly inclined 
smooth surface and driven by a rapidly vibrating bottom 
wall, like in Ref. ^lj. After the "granular gas" reaches a 
steady state, one stops the driving and follows the cool- 
ing dynamics with a fast camera and a particle tracking 
software. While particle rotation and rolling friction may 
prove important, we expect that the main predictions of 
the theory including the scaling behavior of the total 
energy at t — ► t c , will persist. 
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